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Abstract 

A quasi-static (low frequency) model is developed for THUNDER actuators configured 
as displacement sensors based on a simple Raleigh-Ritz technique. This model is used to 
calculate charge as a function of displacement. Using this and the calculated capacitance, 
voltage vs. displacement and voltage vs. electrical load curves are generated and 
compared with measurements. It is shown this model gives acceptable results and is 
useful for determining rough estimates of sensor output for various loads, laminate 
configurations and thicknesses. 
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Introduction 


Thunder actuators are devices constructed from an isotropic laminate of 
aluminum, LaRC-SI adhesive, an electroded Lead Zirconate Titanate (PZT) piezoelectric 
ceramic wafer and a metal backing of either steel or brass. The top aluminum layer is 
used to increase the durability of the laminate and is sometimes omitted. The LaRC-SI 
adhesive is a 1-2 mil thick solid thermoplastic film, and the aluminum layer is indented to 
provide electrical conductivity to the wafer after the actuator is pressure bonded above 
250°C. This hot-melt bonding process induces a mechanical prestress in the THUNDER 
laminate that results a bow shape actuator. This curve is the result of the mismatch in the 
Coefficient of Thermal Expansion (CTE), of the PZT and the metallic prestress layer, as 
the bonded laminate cools to room temperature. If simply supported or arranged in a 
clamshell arrangement the THUNDER wafers produce a linear motion proportional to the 
externally applied voltage. These actuators may also be used as sensors due to the 
intrinsic piezoelectric effect. Since the construction of the THUNDER wafers can be 
varied to a great extent by the selection of materials, processing conditions and geometry, 
it became apparent that computational and analytical models were needed to guide the 
construction of THUNDER wafers to provide practical solutions to engineering 
problems. 

Analytic THUNDER models span more than a decade. The very first was 
developed by Alan R. D. Curtis in December, 1997 who was hired as a NASA contractor 
through BBN Technologies!!]. The BBN model was a very simple linear model that 
assumed only 3 layers, but was a good first effort. In that same timeframe a modeling 
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group at Langley was setup headed by Bill Hunter who was chosen to coordinate the 
modeling efforts of various individuals involved in modeling THUNDER. This group 
consisted of Joel Campbell, Barmac Taleghani [2], and Jill Marlow. Taleghani and 
Marlow were working on finite element models based on NASTRAN and Campbell was 
working on an analytic model. 

In December, 1997 Hunter introduced Campbell to M. W. Hyer's work in 
composite laminates [3-5]. Campbell reformulated Hyer’s work for isotropic media and 
piezoelectric expansion, and thus was able to apply Hyer’s composite modeling to 
THUNDER. Piezoelectric actuation is analogous to thermal expansion so modifying 
those models for piezoelectric layers was not difficult. By January, 1998 Campbell had a 
fully functional model and by February 1998 he submitted two papers to Smart Materials 
and Structures for publication, which were accepted in July on condition of submitting a 
minor revision. The first of these papers detailed the basic model for rectangular actuators 
and the second dealt with a method for predicting blocked force. By March he had a 
model for circular THUNDER actuators and he submitted a patent disclosure, LAR- 
1 5827, based on all 3 papers in April, 1998. He was asked by the NASA patent office not 
to publish those results until an intellectual property judgment could be made, so the 
revision was withheld. Due to a serious backlog this process took over 5 years and by then 
the work was no longer current. Unfortunately, the decision was made not to pursue a 
patent, but disclosed the content to Virginia Power, now Dominion Energy, under their 
THUNDER license agreement. These 3 papers have since been republished as NASA 
TM's [6-8] for historical purposes. 

During this timeframe Hyer was collaborating with a separate NASA group on 
THUNDER. He independently developed an analytical model of his own that he submitted 
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to Smart Materials and Structures in April of 1998 [9], two months after Campbell's 
submission to the same journal. This model was very similar to Campbell’s first paper [6] 
with the following exceptions. First, Hyer formulated his energy integral in terms of forces 
and moments whereas Campbell followed Hyer’s earlier formulation in terms of stresses 
and strains. Although this may look different superficially, the formulation is equivalent. 
Another difference is Campbell’s formulation assumed purely isotropic layers whereas 
Hyer’s model was more general, similar to Hyer's composites work. Campbell also takes 
into account external loading in his model for blocked force. Hyer then followed this up 
with a model for circular actuators in 2002 and 2003 [10-11]. This was a more 
sophisticated model than Campbell’s circular actuator model from 1 998 in that it used 
higher order polynomials. Hyer’s group has since perfected his rectangular actuator model 
to include higher order polynomials and 21 independent coefficients, compared to 4 in 
the earlier model [12], which has been verified experimentally by Mossi et al [ 1 3] . Other 
modeling efforts by others include what appears to be a more sophisticated version of the 
BBN model [14], and some very interesting work modeling hysteresis [15]. 

There has been much less work modeling THUNDER as a sensor, with only two 
published papers using finite element analysis [16-17]. Campbell , however, developed a 
model for THUNDER as a sensor back in 1998 complete with software and examples. 
Unfortunately, due to other commitments, this work was abandoned before the report was 
completely finished. That report was an attempt to model experimental data that had 
been collected by a student a few moths before, and to realize THUNDER actuators as 
potential power sources for energy harvesting. This paper is based on that unpublished 
report. No analytic model for THUNDER actuators as sensors has been published in the 
past decade. This is surprising considering the interest in using THUNDER actuators as 


4 



sensors and energy harvesters. This paper represents a first contribution to that effort, 
which is useful to those who neither have nor wish to use finite element software. 

The Basic Model 

The Von Karman approximation assumes large displacement and small strain. 
This suggests including second order terms in w in addition to the linear terms for u, v 
and w. Expressing the strain to include these terms gives[18], 
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where u 0 , v 0 , and w 0 are the middle thickness displacements in the x, y, and w direction. 
See Figure 1 . 

The total energy takes the form 

C / = jjj {dU 0 - dW T - dW A ) - JJ w p(x,y) dx dy (4) 

where dUo is the stored elastic energy volume density and dWx is the temperature 
contribution. Excluding the pressure term, this follows that used by others [3-5, 19]. The 
last integral is the work done against the external pressure, p. The elastic energy portion 
for the k'th layer can be found with the help of [20] , 
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For the thermal portion we have 
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For the actuation portion we have 
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Actuation is assumed to occur on the j'th layer and Q is the reduced stiffness with 
Ql l=Q22=Ql and Ql2=Q21=Q2- Now that we have the general form for the energy 
integral we may try a Raleigh-Ritz type solution. The method involves making a guess 
for the solution. This guess should look qualitatively like the expected solution yet be 
flexible enough to allow for adjustment. This adjustment comes in the form of changing 
the guess based on a finite number of parameters, which the guess contains. The 
parameters are adjusted in such a way as to minimize the total energy of the system. In 
this situation it would be advantageous to choose a guess that is close enough to the 
classical lamination result that in the linear limit (small scaling) the solution approaches 
the classical lamination result. Such a guess takes the form 
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Classical lamination theory assumes no shearing strain between layers for thermal 
expansion. If we make the same assumptions here we find u° and v° must take the form 
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The ansatz represented by Equation 8 and 9 has been used successfully by others in 
the treatment of laminates and actuators [3-5, 6-8]. Equations 3, 8 and 9 give 
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If the length in the x direction is L x and the length in the y is Ly then, we may 
approximate w 0 ° by, 

„S=-lK, +i ,i4) (ii) 


This assumes the device is simply supported. For p(x,y) we choose a point load of, 


p(x,y) = FS(x)S(y) (12) 


Other choices are possible such as a constant distributed pressure. However, a point 
load is closest to the type of load used with the rectangular devices. The solution is 
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obtained by minimizing the total energy with respect to a, b, c and d. The result of the 
minimi z ation is, 
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The A's, B's and D's have the usual meanings of extensional, coupling and 
bending stiffness. The force and moment per unit width, N' and M', are given by, 
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The z's are measured with respect to the middle surface and zo is the top surface. The 
lengths L x and Ly are the length and width of the actuator. It is a simple matter to show 
that the solution in the linear limit (small scaling - L x =Ly=0) matches the classical 
lamination result if we take a=-K x , b=-Ky, c=e x and d=£y and F=0. 

Unlike the classical lamination result, the non-linear analysis predicts unequal 
curvatures that depend on the magnitude of the scaling and the aspect ratio. It also 
predicts multiple possible solutions. The nature of these solutions are dome-like solutions 
which approach cylinders in the limit of large scaling. Under load, saddles are also 
possible. Typically, there will be one or two stable solutions and possibly a third unstable 
solution depending on the size and material properties. There are also unphysical 
solutions in terms of complex numbers. 
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Figure 1. Experimental setup and Axis definitions. 

As a simple example we take 3 layer actuator constructed from 1 mil brass, 1 mil 
LaRC SI thermoplastic and 6.8 mil PZT 5A. By solving Equation 13 with V=0 and 
plotting F as a function of dome height we find the spring force for the actuator. The 
results of this calculation show the spring force is multivalued. When initially in the short 
axis mode, the actuator may be depressed until a critical load is reached. At this point, the 
actuator switches and locks to the long axis solution. This has been verified 
experimentally and is a general characteristic of all Thunder actuators that have unequal x 
and y curvatures. When the x and y curvatures are strongly unequal, the spring force is 
almost linear except for the bifurcation point. 
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Figure 2. Calculation of spring force vs. distance for a 3-layer actuator. 

Charge and Voltage Generation 

If the device is used as a sensor a charge will be developed across the plates of the 
device. The magnitude of the charge will depend on the stress state of the device. The 
main contributions to this will be from the planar stress due to the mechanical advantage. 
The charge developed from a small piece of the device of dimensions dx dy is, 


With this formulation the strains are functions of x and y and the PZT is assumed to be on 
the j'th layer. By Equation 10 we have, 
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where the strains we calculated are at the center of the PZT layer. The total charge 
generated is obtained by integrating over the surface. The result is, 
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The voltage is obtained by dividing the total charge by the capacitance. The capacitance 
of the device is, 

L L 
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where p is the permittivity of PZT. The voltage is then, 
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To find the voltage one may solve Equation 13 for a, b, c and d in terms of the applied 
force and material properties. It should be noted we have ignored prestress effects on the 
planar charge constant, d 31 . This charge constant is a function of stress as well. 

Voltage and power under load 
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The equivalent circuit for a piezoelectric sensor for low frequencies is a voltage source in 
series with a capacitor. The magnitude of the voltage depends on how the material is 
stressed. For the case of Thunder, this voltage is given by Equation 20. The capacitance is 
the capacitance of the device and is given by Equation 19. 



Figure 3. Equivalent circuit for a Thunder sensor (low frequency model). 

To find the voltage and power under load we first add up all voltages around the loop for 
an external load, R. The result is, 
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Using I=dq/dt and rearranging Equation 21 we find, 
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The steady state solution (solution for long times) of this is, 
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The current is then, 
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The Power may be calculated from, 
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The average power is, 
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By maximizing Equation 26 with respect to R, it can be shown that the power is 
maximized when, 
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This is the point where there is a perfect impedance matching between the Thunder 
sensor and the electrical load, R. 

Experimental 

Voltage vs. displacement: 

To compare this with actual manufactured THUNDERs, a piece measuring 3" x 1 .75" 
was selected and cycled in displacement from 0 to .3 inches at 1 to 10 Hz. The device 
was constructed from 6.8 mil PZT5A, 1 mil LaRC SI and 2 mil brass, and had an initial 
dome height of approximately 0.7 inches. The material properties of the PZT were 
assumed to be d31=-7.83xl0 10 C/lb, Young's modulus of Y=9.1xl0 6 psi, a P=.31 
Poisson's ratio, and a 1700 dielectric constant. The Brass assumed to be Y=1 .5xl0 7 psi 
and P=.34 . The LaRC SI adhesive layer was assumed to be Y=5.8xl0 5 psi and P=.45 . 
The voltage output of the device was measured using an oscilloscope and recorded 
manually. Two samples are shown in the graphs. Figure 4 is a comparison between the 
theoretical curve and the measured data at 1 Hz and 10 Hz. The theoretical curve was 
calculated based only on material properties, geometry and displacements. The 1 Hz 
curve starts to deviate from the ideal curve almost immediately whereas the 10 Hz curve 
matches more closely. The deviation in the 1 Hz curve could be due to creep in the PZT. 
Another possibility would be due to the impedance load of the measuring device, which 
unfortunately is not known due to the fact the data was collected by a summer student in 
1997 who is now long gone. The THUNDER impedance is l/(to C), where C=2.87xl0' 7 


16 



Q . With Z=5.55xl0 3 Q at 1 Hz and Z=5.55xl0 4 Q atlO Hz, a probe in the 1 MQ -20 MQ 


range, could easily cause the effects shown. 



measured (10 Hz) 


Figure 4. Voltage vs. displacement (two samples) 

Voltage, power, and phase vs. electrical load: 

To find power vs. electrical load the sample was cycled in displacement from 0 to 
.15 inches under various resistive loads. The voltage was measured across the resistor and 

the average power was calculated from P=l/2 V^/R. This was then compared to the 
model based only on material properties, geometry and displacements. Here we see the 
deviation increases with increasing resistance. This could either be a material effect or 
possibly due to the impedance of the measuring device. In any even it shows the validity 
of Equation 26. The qualitative nature of the curves is virtually identical and the 
maximum power is obtained by matching the impedance of the sensor with the resistive 
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load. The phase difference between the predicted and measured values indicates the 
impedance of the scope used to measure the voltage is significant. 



Figure 5. Voltage vs. resistance (two samples) 
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Figure 6. Power vs. resistance (two samples). 
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Figure 7. Relative phase vs. resistance (two samples). 


Conclusions 

These results show the basic assumptions used to predict the output of the devices 
are correct and it represents a very simple but useable model for low frequencies. It does, 
in fact, represent the first and only model not based on finite element analysis for the 
prediction of voltage outputs of THUNDER actuators as a function of displacement. 

The impedance of the scope used to measure the voltage output is significant 
enough to affect the voltage output in situations where the load resistance was high, 
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which one would expect. The error in that case was very small for low resistance and 
increased to about 20% over the range measured. A charge amplifier [21] would have 
been better to use in this type of situation since it eliminates the effects of the impedance 
of the scope. The maximum power measurements varied from the model by about 14- 
20% depending on the frequency. Although it is hard to see from the graphs, the 1 Hz 
case was slightly worse that the 10 Hz case. Likewise, the relative phase measurements 
for the 10 Hz case were better than the 1 Hz case as well. 

More work needs to be done to both increase the sophistication of the model and 
to take further measurements using a charge amplifier or high impedance device to 
eliminate the effect of the impedance of the measuring device. 
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